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Abstract 

Many real-world networks, including social and information networks, 
are dynamic structures that evolve over time. Such dynamic networks are 
typically visualized using a sequence of static graph layouts. In addition 
to providing a visual representation of the network topology at each time 
step, the sequence should preserve the mental map between layouts of con- 
secutive time steps to allow a human to interpret the temporal evolution 
of the network. In this paper, we propose a framework for dynamic net- 
work visualization using regularized graph layouts. Regularization encour- 
ages stability of the layouts over time, thus preserving the mental map. The 
proposed framework involves optimizing a modified cost function that aug- 
ments the cost function of a static graph layout algorithm with a grouping 
penalty, which encourages nodes to stay close to other nodes belonging to 
the same group, and a temporal penalty, which encourages smooth move- 
ments of the nodes over time. We introduce two dynamic layout algorithms 
under this framework, namely dynamic multidimensional scaling (DMDS) 
and dynamic graph Laplacian layout (DGLL), that are regularized versions 
of their static counterparts. We apply the proposed algorithms on several data 
sets to illustrate the importance of regularization for producing interpretable 
visualizations of dynamic networks. 



1 Introduction 

The study of networks has emerged as a topic of great interest in recent years, 
with applications in the social, computer, and life sciences, among others. Dy- 
namic networks are of particular interest because networks observed in the real 
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world often evolve over time due to the creation of ne\y node s and edges and 
the re moval of old nodes and edges ( Kossinets and Wattsl 2006 : Leskovec et al. 



20071) . Many developments have been made in mining dynamic ne t works , includ 



ing finding low-rank approximations (ISun et all 120071 : iTong et all l2008h and the 



detec t ion of clusters or communities and how they evolve over time (IChi et al. 



20091 : iMucha et all l2010l : IXu et all. 1201 lal) . However, the closely related task of 



visualizing dynamic networ ks has remained an open problem that has attracted 
attent ion from sociologists (*Mood v et al. , 2005 : Bender-deMoll and McFarland , 
200ddLevdesdorff and Schank, 2008^ 

and the information visualization commu- 



nity ( Brandes and Wagneii 1 19971: iBrandes and Cormanl l2003l : Erten et all l2004l : 
Brandes et al. . 2006; Frishman and Tai 2008 ) among others. Visualization is an 



important tool that can provide insights and intuition about networks that cannot 
be conveyed by summary statistics alone. 

Visualizing static networks is a challenge in itself. Static networks are typically 
represented by graphs, which have no natural representation in a Euclidean space. 
Many graph layout algorithms have been developed to create aesthetically p leasing 

2-D representations of graphs ( Di Battista et al. , 19991 : Herman et al. . 2000). Com- 

mon layout methods for general graphs i nclude force-directed layout (IKamada and KawaiL 



1989';'Fruchterman and Rei ngoldll991l). multidimens ional scaling (MPS) (Ide Leeuw and Heiser , 
1980 : iGansner et al.. .2004. : iBorg and Groenenl. 120050 an d graph Laplacian layout 

' " Koren ll20( 



12003). 



(GLL), also known as spectral layout (iHalll 119701 : 

Dynamic networks are typically represented by a time-indexed sequence of 
graph snapshots; thus visualizing dy namic network s in 2-D presents an additional 
challenge due to the temporal aspect (|BrankeLl2001i) . If one axis is used to represent 
time, then only a single axis remains to convey the topology of the network. While 
it is possible to identify certain trends from a 1-D tim e plot created in this manner , 
it is a poor representation of the network topology. Bran des and Cormanl (120031) 
presented a possible solution to this problem by creating a pseudo-3-D visualiza- 
tion that treats 2-D layouts of each snapshot as layers in a stack. Unfortunately, the 
resulting visualization is difficult to interpret. The more conventional approach is 
to pre sent an animated 2-D layout that evolves over time to reflect the curren t snap- 
shot (Erten et al.', '2004'; 'Moodv et al.', '2005'; 'B ender-deMoll and McFarlandl l2006l : 
Ifirandes et al.. 2006; Frishman^ and TaL 2008 ). A challenge with this approach is to 
preserve the mental map jMisue et al. , 1995 ) between graph snapshots so that the 
transition between frames in the animation can be easily interpreted by a human 
viewer. In particular, it is undesirable for a large number of nodes to drastically 
change positions between frames, which may cause the viewer to lose reference of 
the previous layout. 

Some of the early work on dynamic n etwork visualization simply involved cre- 
ating interpolated transition layouts ((Moody et al.ll2005l : lBender-deMoll and McFarland . 
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20061) . While interpolation does make an animation more aesthetically pleasing, 
it does not constrain the successive layouts in any way to make them more in- 
terpretable. In many real networks, individual snapshots have high variance, so 
creating a layout for each snapshot using a static graph layout method could result 
in large node movements between time steps. Often, this is not due to a failure 
of the static graph layout algorithm but simply a consequence of the cost func- 
tion it is attempting to optimize, which does not consider any other snapshots. 
When dealing with dynamic networks, better performance can be obtained by us- 
ing regularized methods that consider the historical snapshots in addition to the 
current snapshot. Such an approach has been used to develop regularized cluster- 



. 20091: 


Mucha et al. , 


2010; 


Xuetal., 



traditional static clustering algorithms in the dynamic setting. In the context of 
dynamic network visualization, regularization encourages layouts to be stable over 
time, thus preserving the mental map between snapshots. The concept of regular- 
ization has been employed in many problems in statistics and machine learning, 
including ridge regression (Hoerl and Kennard, 1970), also known as Tikhonov 
regularization, the L ASSO ([TibshiraniL 1 1996h . and penahzed matrix decomposition 
dWitten et all l2009h . It is often used to introduce additional information or con- 



straints and to give preference to solutions with certain desirable properties such 
as sparsity, smoothness, or in this paper, dynamic stability in order to preserve the 
mental map. 

We introduce a framework for dynamic network visualization using regularized 
graph layouts. The framework is designed to generate layouts in the on-line setting 
where only present and past snapshots are available, although it can also be used 
in the off-line setting where access to future snapshots is also available. It involves 
optimizing a modified cost function that augments the cost function of a static 
graph layout algorithm with two penalties: 

1. A grouping penalty, which discourages nodes from deviating too far from 
other nodes belonging to the same group; 

2. A temporal penalty, which discourages nodes from deviating too far from 
their previous positions. 

Groups could correspond to a priori knowledge, such as participant affiliations in 
social networks. If no a priori group knowledge is available, groups can be learned 
from the network using, for example, the aforementioned evolutionary clustering 
algorithms. The grouping penalty keeps group members close together in the se- 
quence of layouts, which helps to preserve the mental map because nodes belong- 
ing to the same group often evolve over time in a similar fashion. The tempo- 
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ral penalty helps to preserve the mental map by discouraging large node move- 
ments that may cause a human to lose reference of previous node positions, such 
as multiple nodes moving unnecessarily from one side of the layout to the opposite 
side. Using the proposed framework, we develop two dynamic layout algorithms, 
dynamic multidimensional scaling (DMDS) and dynamic graph Laplacian layout 
(DGLL), that are regularized versions of MDS and GLL, respectively. 

The idea of preserving temporal stability in dynamic graph layouts has been 

proposed in previous wo r k (IBrandes and Wagneiill997l : lErten et al.l.l2004l : lBaur and Schankl . 



20081 : iFrishman and Tai |2008|). The idea of placing nodes b elonging to the same 
group together in a layout has also appeared in previous w ork (IWang and Miyamotol 



19961 : iHuang and Eadesl. 119981 : ICosta and Hero 1111120051 ). albeit in the static graph 



setting. To the best of our knowledge, this is the first work that presents a frame- 
work for dynamic network visualizatior0 that incorporates both grouping and tem- 
poral regularization. The algorithms DMDS and DGLL are also unique in that 
they add both grouping and temporal regularization to M PS and GLL, respe c tively . 
Existing methods for temporal regul arization in MDS (Baur and Schank . 2008h 
and grouping regularization in GLL (ICosta and Hero IIll |2005) are subsumed by 
DMDS and DGLL, respectively. The methods for grouping regularization in MDS 
and temporal regularization in GLL used in this paper are novel. We apply the pro- 
posed algorithms on a selection of dynamic network data sets to demonstrate the 
importance of both grouping and temporal regularization in creating interpretable 
visualizations. 



2 Background 

We begin by specifying the notation we shall use in this paper. Time-indexed 
quantities are indicated using square brackets, e.g. X[t]. We represent a dynamic 
network by a discrete-time sequence of graph snapshots. Each snapshot is repre- 
sented by a graph adjacency matrix W[t] where Wij[t] denotes the weight of the 
edge between nodes i and j at time t (chosen to be 1 for unweighted graphs), and 
Wij[t] = if no edge is present. We assume all graphs are undirected, so that 
Wij[t] = Wji[t]. For simplicity of notation, we typically drop the time index for all 
quantities at time step t, i.e. W is assumed to denote W[t]. 

We refer to a graph layout by a matrix X € M"^^, where n denotes the number 
of nodes present at time t, and each row X(j) corresponds to the s-dimensional 
position of node i. We are most interested in 2-D visualization {s = 2), although 
the proposed methods can also be applied to other values of s. The ith column of 
X is denoted by Xj, and the individual entries by Xij. The superscript in xj'^^ is 

A preliminary version of this work can be found in jXu et alll20TTbh . 
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used to denote the value of Xj at iteration h of an algorithm. The norm operator 
II • II refers to the /2-norm, and tr(-) denotes the matrix trace operator. We denote 
the all ones column vector by 1. 

We now summarize the static graph layout methods of multidimensional scal- 
ing and graph Laplacian layout, which operate on a single graph snapshot. We 
develop regularized versions of these methods for dynamic networks in Section [3] 



2.1 Multidimensional scaling 

Multidimensional scaling (MDS) is a family of statistical methods that aim to find 
an optimal layout X G M"^^ such that the distance between X(j) and xq) for all 
i / j is as close as possible to a desired distance d-ij. There are a variety of different 
cost functions and associ ated algorithms for MDS; we refer interested readers to 



Borg and GroenenI (120051) . Here we describe the cost function known as stress and 



its associated majorization algorithm. The stress of a layout X is given by 

Stress(X) = - ^ {6ij - - X(j)||)^ , (1) 

i=i j=i 

where Vij > denotes the weight or importance of maintaining the desired distance 
6ij. We refer to the matrix V = [vij] as the MDS weight matrix to avoid confusion 
with the graph adjacency matrix W, which is also sometimes referred to as a weight 
matrix. In order to use stress MDS for graph layout, the graph adjacency matrix W 
is first converted into a desired distance matrix A = [Sij]. This is done by defin- 
ing a distance metric over the graph and calculating distances between all pairs of 
nodes. The distance between two no des i and 7 is typically taken to be the length of 
the shortest path between the nodes ( Gansner et al. , 200 4^. For weighted graphs, it 



is assumed that the edge weights denote dissimilarities; if the edge weights instead 
denote similarities, they must first be converted into dissimilarities before comput- 
ing A. The MDS weights Vij play a crucial role in the aesthet ics of the layout. The 



commonly used Kamada-Kawai (KK) force-directed layout (IKamada and KawaiL 



19891 ') is a special case of stress MDS where the weights are chosen to be Vij = 5^^ 



for i / i and vu = 0. 

The objective of stress MDS is to find a layout X that minimizes ([T]). ^ can 
be decomposed into 



^ n n ^ n n n n 

i=l j=l i=l i=l i=l i=l 

Note that the first term of © is independent of X. The second term of ^ can be 
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written as tT{X'^ RX) where the n x n matrix R is given by 



i = j. 



(3) 



tT{X^ RX) is quadratic and convex in X and is easily optimized. 

The third term of (ID cannot be written as a quadratic form. However, it can be 
optimized by an iterative majorizat ion method known as "scaling by majorizing a 
complicated function" (SMACOF) ( de Leeuw and Heiser . 198Cll : Borg and Groenen , 
This non-quadratic term is iteratively majorized, and the resulting upper 
bound for the stress is then optimized by differentiation. For a matrix Z E M"^*, 
define the matrix-valued function S{Z) by 



Sij{Z) 



- Z(i)ll i^j, 

Then, it is shown in (iGansner et al. I I2OO4I : iBorg and GroenenL 120051) that 



(4) 



n n 

tr{X^S{Z)Z)>--J2Y. 



^ij^ij ll-^(j) 



i=i j=i 



so that an upper bound for the stress is 



2 E E + tr(X^i?X) - 2 tr(X^5(Z)Z). 

i=i j=i 



(5) 



By setting the derivative of Q with respect to X to 0, the minimizer of the upper 
bound is found to be the solution to the equation RX = S{Z)Z. 

The algorithm for optimizing stress is iterative. Let X'^'^^ denote an initial lay- 
out. Then at each iteration h, solve 



(6) 



Ah) 



for Xa"'' for each a = 1, . . . , s. Q can be solved using a standard linear equation 
solver. Note that R is rank -deficient; this is a co nsequence of the stress function 
being translation-invariant (IGansner et all |2004[) . The translation-invariance can 
be removed by fixing the location of one point, e.g. by setting X(i) = 0, removing 
the first row and colu mn of R, and removing the first row of S'(X(^^^))X(''^^) 
(IGansner et al.ll2004b . (|6]l can then be solved efficiently using Cholesky factoriza- 
tion. The iteration can be terminated when 

Stress (XC^-I)) - Stress {X'^^y) 



< e, 



Stress (XC^-i)) 
where e is the convergence tolerance of the iterative process 



(V) 
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2.2 Graph Laplacian layout 



Graph Laplacian layout (GLL) methods optim ize a quadratic function associated 
with the graph Laplacian matrix (lKorenLl2003l) . which we call the GLL energy. The 
graph Laplacian is obtained from the adjacency matrix hy L = D — W, where D 
is the diagonal matrix of node degrees defined by da = Yl^=i "^ij- ^'^^ weighted 
graphs, GLL assumes that the weights correspond to similarities between nodes, 
rather than dissimilarities as in MDS. GLL is also referred to as "spectral layout" 
because the optimal solution involves the eigenvectors of the Laplacian, as we will 
show. The GLL energy function is defined by 



^ n n 

Energy(X) = i^^^ mj ||x(i) - x 



(8) 



It is easily shown that Energy (X) = tr(X^LX). The GLL problem can be for- 
mulated as (Eaii 1 1 970l : TKorenL boOBh : 



mill ti(X^ LX) 

X 

subject to X^X = nl 
X'^1 = 0. 



(9) 
(10) 
(11) 



From dH]), it can be seen that minimizing the GLL energy function aims to make 
edge lengths short by placing nodes connected by heavy edges close together. ([TTI ) 
removes the trivial solution = 1, which places all nodes at the same location in 
one dimension. It can also be viewed as removing a degr ee of freedom in the layout 
due to translation invariance (|Belkin and NiyogiL |2003|) by setting the mean of x^ 
to for all a. Since x^ has zero-mean, (x^Xa)/n corresponds to the variance 
or scatter of the layout in dimension a. Thus (flOl ) constrains the layout to have 
unit variance in each dimension and zero covariance between dimensions so that 
each dimension of the layout provides as much additional information as possible. 
M oreover, one can see that (flOl ) diffe rs sU ghtly from the usual constraint X^ X = 
/ ("Be lkin and NiyogiL I2OO3. : Xoreiil 12003 ). which constrains the layout to have 
variance 1 /n in each dimension. In the dynamic network setting where n can vary 
over time, this is undesirable because the layout would change scale between time 
steps if the number of nodes changes. 



By a generalization of the Rayleigh-Ritz theorem (ILutkepohll . 119971) . an opti- 
mal solution to the GLL problem is given by X* = -v/n[v2, V3, . . . , v^+i], where 
Vj denotes the eigenvector corresponding to the ith smallest eigenvalue of L. Note 
that vi = {l/^/n) \ is excluded because it violates the zero-mean constraint (ITTI ). 
Using the property that iT{ABC) = iT{CAB), the cost function Q is easily 
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shown to be invariant to rotation and reflection, so X*Ris also an optimal solution 
for any R^R = I. 

In practice, it has been found that using degree- normalized eigenv ectors of- 
ten re sults in more aesthetically pleasing layouts ( BeMn and Niyogi , 2003. ; Korenj 



20031) . The degree-normalized layout problem differs only in that the dot product in 
each of the constraints is replaced with the degree-weighted dot product, resulting 
in the following optimization problem: 



min triX^ LX) 
X 



subject to triX'^DX) = tr{D)I 
X'^Dl = 0. 



The optimal solution is given by X* = yM^rp?) [u2, U3, . . . , u^+i] or any rotation 
or reflection of X*, where Uj denotes the generalized eigenvector corresponding to 
the ith smallest generalized eigenvalue of {L,D). Again, ui = (l/y^ ii{D)) 1 is 
excluded because it violates the zero-mean constraint. A discussion on the merits 
of the degree normalization can be found in ( Korenj 2003). 



3 Regularized layout methods 
3.1 Regularization framework 

The aforementioned static layout methods can be applied snapshot-by-snapshot to 
create a visualization of a dynamic network; however, the resulting visualization 
is often difficult to interpret, especially if there are large node movements between 
time steps. We propose a regularized layout framework that uses a modified cost 
function, defined by 

C-modified — Cgtatic ~l~ CtCgj-ouping ~l" /^Ctemporal- (12) 

The static cost Cstatic corresponds to the cost function optimized by the static layout 
algorithm. For example, in MDS, it is the stress function defined in ([T]l, and in GLL, 
it is the energy defined in The grouping cost Cgiouping is chosen to discourage 
nodes from deviating too far from other group members; a controls the importance 
of the grouping cost, so we refer to aCgrouping as the grouping penalty. Similarly, 
the temporal cost Ctemporai is chosen to discourage nodes from deviating too far 
from their previous positions; /3 controls the importance of the temporal cost, so 
we refer to /3Ctemporai as the temporal penalty. We propose quadr atic forms for these 
penalties, similar to ridge regression dHoerl and Kennard , 1970 ). 
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Let k denote the number of groups. Define the group membership by an n x A; 
matrix C where 




1 node z is in group I at time step t, 
otherwise. 



We introduce grouping regularization by adding group representatives, which also 
get mapped to an s-dimensional position, stored in the matrix Y G R*^^*. The 
proposed grouping cost is given by 

k n 

Cgroupmg(^,^) =^^Cii\\x^i^-yQ)f, (13) 

1=1 i=l 

where y(;^ denotes the position of the Ith representative. Notice that the grouping 
cost is decreased by moving y(;) and X(j) towards each other if node i is in group I. 
As a result, nodes belonging to the same group will be placed closer together than 
in a layout without grouping regularization. Notice also that we do not require 
knowledge of the group membership of every node. Nodes with unknown group 
memberships correspond to all-zero rows in C and are not subject to any grouping 
penalty. 

We introduce temporal regularization on nodes present at both time steps t and 
i — 1 by discouraging node positions from deviating significantly from their pre- 
vious positions. This idea is often referred to in the information visualization ht- 
erature as maintaining dynamic stability of the layouts and is often used to achieve 
the goal of preserving the mental map. Define the diagonal matrix E by 




1 node i was present at time step t — 1, 
otherwise. 



The proposed temporal cost is then given by 

n 

Ctemporal(^,^[i - 1]) = ^eii||x(,) - X(i)[t - (14) 

i=l 

The temporal cost is decreased by moving X(j) towards X(j) [t—1], but unlike in the 
grouping cost, X(j) [t — 1] is fixed because it was assigned at the previous time step. 
Thus the previous node position acts as an anchor for the current node position. 

We note that the temporal cost measures only the stability of the layouts over 
time and is not necessarily a measure of goodness-of-fit with regard to the dynamic 
network. For example, if there is a sudden change in the network topology, an 
extremely low temporal cost may be undesirable because it could prevent the layout 
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from adequately adapting to reflect this change. Thus one must consider the trade- 
off of adaptation rate versus stability when choosing the temporal regularization 
parameter. 

Next we demonstrate how the grouping and temporal penalties can be intro- 
duced into MDS and GLL as examples of the proposed regularization framework. 



3.2 Dynamic multidimensional scaling 

The dynamic multidimensional scaling (DMDS) modified cost is given by the mod- 
ified stress function 



MStress(X,y) = - '^^Vij{5ij - ||x(i) -^{^j)\\f 
i=i j=i 

k n n 

+ a X] X] - y(fc) f + /3 ^ ||x(i) - X(j) [t 



(15) 



1=1 i=l 



i=l 



The first term of (fTSl ) is the usual MDS stress function, while the second term 
corresponds to the grouping penalty, and the third term corresponds to the tempo- 
ral penalty. The constants a and /3 are the grouping and temporal regularization 
parameters, respectively. 

To optimize ([T5] ). we begin by re-writing the first two terms into a single term. 
Define the augmented MDS weight matrix by 



V 



V aC 
aC^ 



(16) 



where the zero corresponds to an appropriately sized all-zero matrix. Similarly, 
define the (n + A;) x (n + k) augmented desired distance matrix A by filling the 
added rows and columns with zeros, i.e. 



Let 



X 



A 0' 


X 
Y 



(17) 



(18) 



denote the positions of the both the nodes and the group representatives. Then, the 
first two terms of ([T5] ) can be written as 



n+fc n+fc 



i=i i=i 



Hi) 



X(4 
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which has the same form as the usual stress defined in ([T]). The third term in ([Ts] ) 
can be written as a quadratic function of X, namely 



/3 trlX^ EX) - 2tr[X'^EX[t - 1]) + tr (x'^[t - l]EX[t - 1] 



where the (n + A;) x (n + k) matrix E and the {n + k) x s matrix X[t — 1] are 
constructed by zero-filling as in the definition of A. 

Following the derivation in Section ITTl for any {n + k) x s matrix Z, ( fT5] ) can 
be majorized by 



n+k n+k 

2 

1=1 j=i 



- V V Vifdfj + iT{X^RX) -2tT{X^S{Z)Z) + p\tT{X^EX) 
^=l i=i ^ (19) 



- 2 tr(X^ EX[t-l])+ tT{X^ [t - l]EX[t - 1]) 

where R and 5 are defined by substituting the augmented matrices V and A for V 
and A, respectively, in (O and ( [T9l ) is quadratic and convex in X so the min- 
imizer of the upper bound is found by setting the derivative of (fT9l) to 0, resulting 
in the equation 

{R + I3E)X = S{Z)Z + pEX[t - 1]. 

This can again be solved sequentially over each dimension. As in Section 12.11 we 
solve this iteratively using the previous iteration as the majorizer, i.e. at iteration 
h, solve 

{R + f3E)5S^^ = S (x('^-i)) 4'^-!) + pES^a[t - 1]. (20) 



for xi'^'' for each a = 1, . . . , s. The process is iterated until the convergence crite- 
rion ^ is attained. The first iterate can be taken to be simply the previous layout 
Xa[t — 1]. Unlike in ordinary MDS, the system of linear equations in (l20l ) has a 
unique solution provided that at least a single node was present at time step t — 1, 
because R + (5E has full rank in this case. 

Pseudocode for the DMDS algorithm for t = 1,2,... is shown in Fig.[T] where 
shortest_paths(-) computes the matrix of shortest paths between all pairs of nodes, 
and MDS_weights(-) computes the MDS weight matrix. (l20l ) can be solved by 
performing a Cholesky factorization on {R + fiE) followed by back substitution. 
At the initial time step (t = 0), there are no previous node positions to initialize 
with, so a random initialization is used. Also, the position of one node should be 
fixed before solving (l20l) due to the translation-invariance discussed in Section [ZT] 
The time complexity of the algorithm at all subsequent time steps is dominated 
by the O(n^) complexity of the Cholesky factorization, assuming k ^ n, but 
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1: fort = 1,2,... do 

2: A ^ shortest_paths(M^) 

3: V ^ MDS_weights(A) 

4: Construct V and A using ([T6l ) and ([TTl) . respectively 

5: Construct R by substituting y for F in ([3]) 

6: /l ^ 

7: ^ X[t - 1] 

8: repeat 

9: h + 1 

10: Construct S'(x(^-i)) by substituting y, A, and xC'"^) for A, and Z, 

respectively, in ^ 
11: for a = 1, . . . , s do 

12: Solve {R + l5E)^f^ = ^(x('*-i))xi''"^^ + l3E^a[t - 1] for ^^a^ 

13: end for 

14: until [ MStress (X^^-^)) - MStress (X^'^))] / MStress (X^'^-i)) < e 

15: X ^ X^'^) 

16: end for 

Figure 1 : Pseudocode for the DMDS algorithm. 

the factorization only needs to be computed at the initial iteration {h = 1). All 
subsequent iterations require only matrix- vector products and back substitution and 
thus have O(n^) complexity. 

3.3 Dynamic graph Laplacian layout 

The dynamic graph Laplacian layout (DGLL) modified cost is given by the modi- 
fied energy function 



^ n n 

MEnergy(X,y) = - ^^z/;y||x(,) - X(j)f 



2 

i=l j=l 

(21) 

k n n 
1=1 i=l i=l 

Like with DMDS, the first term of (|2TI ) is the usual GLL energy function, while the 
second term corresponds to the grouping penalty, and the third term corresponds 
to the temporal penalty. Again, the parameters a and /3 correspond to the grouping 
and temporal regularization parameters, respectively. 
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We first re-write (|2TI ) in a more compact form using the grapli Laplacian. De- 
fine tlie augmented adjacency matrix by 



W 



■ W aC 
aC^ 



(22) 



Notice that the group representatives have been added as nodes to the graph, with 
edges between each node and its associated representative of weight a. Define 
the augmented degree matrix by D by da = Yl^=i '^ij-' ^^i^ the augmented graph 
Laplacian by L = Z) — W. The first two terms of (|2TI ) can thus be written as 
ii{X'^LX), where X is as defined in (1181) . The third term of (1211) can be written 
as 



/3 



iT{X^ EX) - 2tT{X'^EX[t - 1]) + iT{X^\t - l]EX[t - 1]) 



(23) 



where E is zero-filled as described in Section 13.21 The final term in (1231) is inde- 
pendent of X and is henceforth dropped from the modified cost. 

We now consider the constraints, which differ depending on whether the lay- 
out is degree-normalized, as discussed in Section 12.21 We derive the constraints 
for the degree-normalized layout; the equivalent constraints for the unnormalized 
layout can simply be obtained by replacing D with the identity matrix in the deriva- 
tion. First we note that, due to the temporal regularization, the optimal layout is no 
longer translation-invariant, so we can remove the zero-mean constraint. As a re- 
sult, the variance and orthogonality constraints become more complicated because 
we need to subtract the mean. Denote the degree-weighted mean in dimension a 
by 



n+k 



fJ-a 



^n+k 



-1 1 



Z^j=l i=i 

Then the degree-weighted covariance between the ath and 6th dimensions is given 
by 

^ n+k 

COv(Xa,Xb) = - ^ dii{xia " fla){Xib - fib) 



n+k 



X^n+k J diiXiaXih ~ \ 

ti{b) ' 



(n+k 
/ , diiXia 
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where M is the centering matrix defined by 



M = D 



tT{D) 



(24) 



Combining the modified cost function with the modified constraints, the nor- 
mahzed DGLL problem is as follows: 



min tr{X^ LX) + f3 tr{X^ EX) - 2X^ EX[t - 1]) 
X L 

subject to tr(X^MX) = tr{D)I. 



(25) 
(26) 



Again, the unnormalized problem can be obtained by replacing D with the identity 
matrix in (l24l ) and ( [26l ). Note that ( |25] ) contains a linear term in X. Hence the 
optimal solution cannot be obtained using scaled generalized eigenvectors as in the 
static GLL problem. (|25l) can be solved using standard algorithms for constrained 
nonlinear optimization (IBazaraa et all 120061) . The cost function and constraints 
consist only of linear and quadratic terms, so the gradient and Hessian are eas- 
ily computed in closed form (see Appendix [Ajl. Unfortunately, the problem is not 
convex due to the equality constraints; thus a good initialization is important. The 
natural choice is to initialize using the previous layout X'^'^'^ = X[t — 1]. To avoid 
getting stuck in poor local minima, one could use multiple restarts with random 
initialization. 

Pseudocode for the DGLL al gorithm for t = 1 , 2, ... is shown in Fig. |2] We 



use the interior-point algorithm of iByrd et all (|1999|) to solve (1251 ). We find in prac- 
tice that random restarts are not necessary unless /3 is extremely small because the 
temporal regularization penalizes solutions that deviate too far from the previous 
layout. For other choices of /3, we find that the interior-point algorithm indeed 
converges to the global minimum when initialized using the previous layout. The 
most time-consuming operation in solving (1251 ) consists of a Cholesky factoriza- 
tion, which must be updated at each iteration. At the initial time step {t = 0), there 
are no previous node positions, and hence, no linear term in (|25] ). so the layout is 
obtained using scaled generalized eigenvectors, as described in Section 12.21 The 
time complexity at all subsequent time steps is dominated by the 0{n^) complexity 
of the Cholesky factorization. 



3.4 Discussion 

We chose to demonstrate the proposed framework with MDS and GLL; however, it 
is also applicable to other graph l ayout methods, such as the Frucht erman-Reingold 
method of force-directed layout ((Fruchterman and ReingoldLll99lh . Since the static 
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1: fort = 1,2,... do 

2: Construct W using (l22l ) and its corresponding Laplacian L = D — W 
3: Construct the centering matrix M using (|24] ) 

4: ^ X[t - 1] 

5: Solve (|25] ) using the forms for V/, (7, and J in Appendix lAl 

6: for r = 1 — )• max -restarts do {if multiple random restarts are necessary} 

7: Randomly assign 

8: Solve (l25l) using the forms for V/, 51, -fT, and J in Appendix lAl 
9: end for 

10: X ^ best solution to ( |25] ) over all initializations 

11: end for 



Figure 2: Pseudocode for the DGLL algorithm. 



cost functions of MDS and GLL encourage different appearances, the same is true 
of DMDS and DGLL. Ultimately, the decision of which type of layout to use de- 
pends on the type of network and user preferences. Kamada-Kawai MDS layouts 
are often preferred in 2-D because they discourage nodes from overlapping due to 
the large MDS weights assigned to maintaining small desired distances. On the 
other hand, if a 1-D layout is desired, so that the entire sequence can be plotted as 
a time series, node overlap is a lesser concern. For such applications, DGLL may 
be a better choice. 

Another decision that needs to be made by the user is the choice of the param- 
eters a and {3, which can be tuned as desired to create a meaningful animation. 
Unlike in supervised learning tasks such as classification, there is no ground truth 
in visualization so the selection of parameters in layout methods is typically done 
in an ad-hoc fashion. Furthermore, multiple layouts created by differing choices 
of parameters could be use ful for visuahzing different p ortions of the network or 
yielding different insights dwitten and Tibshiranil I2OI ih . This is particularly true 
of the grouping regularization parameter a. When a high value of a is used, nodes 
belonging to the same group are placed much closer together than nodes belong- 
ing to different groups. The resulting visualization emphasizes node movements 
between groups (for nodes that change group between time steps) while sacrificing 
the quality of the node movements within groups. On the other hand, when a low 
value of a is used, node movements within groups are more clearly visible, but 
node movements between groups are more difficult to see. We explore the effect of 
changing parameters on the resulting animation in several experiments in Section 

m 
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4 Related work 



The regularized graph layout framework proposed in this paper utilizes a grouping 
penalty that places nodes belonging to the same group together and a temporal 
penalty that places nodes near their positions at neighboring time steps. Node 
grouping and temporal stability in the context of graph layout have previously been 
studied independently of each other. We summarize relevant contributions to both 
of these areas. 



4.1 Node grouping 

Several grouping techniques for static graph layo ut have been previous l y been pro 



posed. Given a partition of a graph into groups, IWang and Miyamoto! (119961) pro- 
pose a modified force-directed layout that considers three types of forces: intra- 
forces, inter-forces, and meta-forces. Intra-forces and inter-forces denote forces 
between nodes in the same group and nodes in different groups, respectively. Meta- 
forces correspond to forces between groups; nodes in the same group are subject to 
identical meta-forces. By decreasing the strength of inter-forces and increasing the 
strength of meta-forces, nodes belonging to the same group get positioned closer 
to each other in the la y out. 



Huang and EadesI (| 19981) developed a system called DA-TU for visualizing 
groups in large static graphs. It also utilizes a modified force-directed layout with 
inter- and intra-forces. However, rather than using meta-forces, DA-TU adds a vir- 
tual node for each group with a virtual force between each virtual node and each 
node in its group. Notice that the virtual nodes are identical to the group represen- 
tatives in our proposed framework; however, the use of virtual forces to achieve 
grouping regularization differs from the squared Euclidean distance grouping cost 
we propose. In addition, DA-TU was designed for visualizing static graphs at many 
scales rather than visualizing dynamic graphs, so it does not contain a temporal sta- 
bility penalty. 

Grouping techniques have been applied in the field of supervised dimension- 
ality reduction, which is very closely related to graph layout. The objective of di- 
mensionality reduction (DR) is to find a mapping (f) -.W , p > s from a high- 
dimensional space to a lower-dimensional space while preserving r nany of the char- 



acteri stics of the data representation in the high-dimensional space (i Lee and Verleysen , 



20071) . For example, MDS is a DR method that attempts to preserve pairwise dis- 
tances between data points. In the supervised DR setting, one also has a priori 
knowledge of the grouping structure of the data. Supervised DR methods pose 
the additional constraint that data points within the same group should be closer 
together in the low-dimensional space than points in separate groups. Notice that 
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this is the same grouping c onstra int we pose in our regularized layout framework. 



Witten and Tibshiranil (1201 ll) proposed a supervised version of MDS (SMDS) 



that optimizes the following cost function over X: 



^ n n 



X(j)||)2 + a 



=1 i=i 



E 

t,j-yj>yi 



a=l 



(27) 

where is an ordinal value denoting the group membership of data point i. Notice 
that the first term in (|27] ) is the ordinary MDS stress with Vij = 1 for all 
while the second term provides the grouping regularization. a controls the trade- 
off between the two terms. The key difference between the SMDS grouping penalty 
and the DMDS grouping penalty proposed in this paper is in the way groups are 
treated. SMDS assumes that groups are labeled with an ordinal value that allows 
them to be ranked, and the form of the grouping penalty in (1271 ) does indeed tend to 
rank groups in R** by encouraging Xja > Xia, a = 1, . . . , s for alH, j : yj > yi. On 
the other hand, our proposed grouping penalty treats group labels as categorical. It 
does not rank groups in W but simply pulls nodes belonging to the group together. 

Another related method for supervised DR is classi fication constrained dimen- 
sionality reduction (CCDR) dCosta and Hero IIll l2005r) . which is a supervised ver- 
sion of Laplacian eigenmaps (|Belkin and NiyogiLl2003h . CCDR optimizes the fol- 
lowing cost function over (X, Y): 



^ n n 



X 



(i)l 



y(/)l 



=1 j=i 



1=1 i=l 



Notice that this cost function is simply the sum of the GLL energy and the DGLL 
grouping penalty. Indeed, DGLL can be viewed as an extension of CCDR to time- 
varying data. The CCDR solution is given by the matrix of generalized eigenvec- 
tors U = [u2, . . . , Us+i] of (L, D), where L and D are as defined in Section [331 
Although the addition of the temporal regularization due to the anchoring presence 
of the previous layout X[f — 1] prevents the DGLL problem from being solved us- 
ing generalized eigenvectors, it discourages large node movements between time 
steps in order to better preserve the mental map. 



4.2 Temporal stability 

There have been many previous studies on the problem o f laying out dynamic net- 
works while preserving stability between time snapshots. iMoody et all (120051) pro- 
posed to generate dynamic layouts using a static layout method such as Kamada- 
Kawai MDS and to initialize at each time step using the layout generated at the 
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previous time step. The idea of this approach is to anchor the nodes initially so that 
the entire layout does not get rotated. The anchoring differs from the temporal reg- 
ularization proposed in this paper, which penalizes changes in node positions over 
time and can be thought of as anchoring all i terations rather than just the initial it- 
eration. The approach of iMoodv et all (120051) is implemented in the s ocial network 
visualization software SoNIA (IBender-deMoU and McFarlandll201 ih . Experimen- 
tally, we find that solely anchoring the initialization is insufficient at preventing 
drastic node movements over time (see Section|5]for examples). 



To enforce stability in layouts at consecutive time steps, iBrandes and Wagner 



(Il997h proposed a probabilistic framework for dynamic network layout where the 
objective is to choose the layout at a particular time step with maximum posterior 
probability given the previous layout. The layout with maximum posterior prob- 
ability can be found by maximizing the product of the prior probability and like- 
lihood of the layout. Similar to our proposed framework, the probabilistic frame- 
work is applicable to a wide class of graph layout methods. The prior probability of 
a layout is taken to be an (appr opriately normalized) invers e exponential function 
of the static cost of the layout. iBrandes and Wagner (Il997h proposed several dif- 
ferent types of likelihood functions corresponding to different notions of temporal 
stability. One such notion is to demand stability of node positions in consecutive 
layouts. For this notion of stability, the authors propose the likelihood for each 
node to be a spherical Gaussian distribution centered at the node's position in the 
layout; the likeUhood of the layout is then taken to be the product of the likelihoods 
of all nodes present in consecutive time steps. Thus the posterior probability can 
be written as (up to a normalizing constant). 



exp 



+ 



El 



2a2 




(28) 



where cr is a scaling parameter for the amplitude of node movements. Notice that 
by taking the logarithm of (|28] ). one obtains the same form as our proposed regular- 
ized framework (fT2] ). excluding the grouping cost, with /3 = 1/(2(7^). The layout 
that maximizes the logarithm of (1281 ) is the same layout that maximizes posterior 
probability; thus, under thi s notion of stability, there is a n equivalence between the 
probabilistic framework of IBrandes and Wagnen (| 19970 and the temporal regular- 
ization framework proposed in this paper. 

Other methods for preserving tem poral stability in dynami c layouts tend to be 
specific to a particular layout method. iBaur and SchankI (l2008h proposed a tempo- 
rally regularized MDS algorithm that uses the following localized update rule at 
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each iteration h for each node i at each time step t: 



(h) _ + /3{eiiXia[t - 1] + fiiXjalt + 1]) 



(29) 



where 



(h-i) _ (h-i) 



^(i) 



and F is the diagonal matrix defined by 

fa ' 



1 node i is present at time step t + 1, 
otherwise. 



This algorith m is implemente d in the soci al netwo rk analysis and vi sualization soft- 

ware Visone ( Brandes and Wa gner. 2 004i : lvisone ) and was used in ( Leydesdorff and Schank . 
2008h for visualizing similarities in journal content over time. ( [29l ) is an off-line 
update because it uses both the node positions at time steps t — 1 and t + 1 to 
compute the node position at time step t, whereas the methods we propose, includ- 
ing DMDS, are on-line methods that use only current and past data. (l29l) can be 
mo dified into an on-line up date by removing the terms involving fa. It was shown 
in ( Baur and Schank . 2008h that the localized update of (l29l ) optimizes the sum of 
the MStress function in (fTSl) over all t with A; = 0, i.e. without a grouping penalty. 
Likewise, the on-line modification optimizes the MStress function at a single time 
step with k = 0. Hence the proposed DMDS layout method ca n be viewed as an 
on-line modification of the method of ( Baur and Schankl . boosl) with the addition 
of a grouping penalty. 

When it comes to GLL, to the best of our knowledge, there is n o prior work that 
explicitly penalizes large node movements. iBrandes et al.l (|2006|) suggest two ap- 
proaches for layout of dynamic networks using GLL. The first is to simply compute 
a static layout for each graph snapshot using the spectral method describ ed in Sec- 
tionl2.2l The eigenvectors are calculated using power iteration ( Trefethen and Bau III , 
19971') initiali z ed wi th the layout at the previous time step, similar to the approach of 
Moody et all ((20051) for MDS. As previously mentioned, this does not ensure lay- 
out stability. The second approach is to compute a spectral layout of a smoothed 
graph Laplacian matrix XL[t — 1] + (1 — X)L\t]. A similar approac h has been used 
for clustering dynamic networks ( Chi et al. , 20091 : Xu et al. , 2011a ). This approach 
incorporates information from two graph snapshots and should perform better than 
the first approach, but it also does not explicitly constrain the layout at time t from 
deviating too far from the layout at time t — l. We henceforth refer to this approach 
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as the BFP method and compare its performance to our proposed DGLL algorithm 
in Section [5] 



Ot her methods for layout of dynamic n etworks have also been proposed (lErten et al 



20041 : iFrishman and Tall l2008h . TGRIP (lErten et all 12004 ) is a modified force 



directed layout method with added edges between vertices present at multiple time 
steps. The user-selected weights of these added e dges control the amount of tem- 
poral regularization in the layouts. The method of lprishman and Tall ( 2008 ) is also 
a modified force-directed layout. It is an on-line method that uses pinning weights 
to previous node positions to achieve temporal regularization and a GPU-based im- 
plementation to reduce run-time. The emphasis in both methods is on improving 
scalability to deal with extremely large networks by coarsening graphs to com- 
pute an initial layout then applying local refinements to improve the quality of the 
layout. As a result, they are applicable to much larger networks than the O(n^) 
methods proposed in this paper. However, these methods do not incorporate any 
sort of grouping regularization to discourage nodes from deviating too far from 
other nodes in the same group. 



5 Experiments 



We demonstrate the proposed framework by applying DMDS and DGLL on a sim- 
ulated data set and two real data sets. Several snapshots of the resulting visual- 
izations are presented. The full, animated visuali zations over all time steps can be 
found on the supporting website dXu et al.ll2012l) . 

In the second experiment, we do not have a priori group knowledge. Hence 
we learn the grou ps using the AFFECT evolutionary spectral clustering algorithm 
(|Xu et al.Ll2011ah . summarized in AppendixE] In the other two experiments, we do 
have a priori group knowledge. We compute layouts both using the known groups 
and the groups learned by clustering. We also compute layouts usin g several ex- 



isting methods for comparison. DM DS is compare d to the method of iMoody et al. 



(120051) used in So NIA_( Bender-deMoll and McFarland. 2011) and the method of 
Baur and SchankI (|2008|) used in Visone ([Brandes and Wagneiil2004l : lvisoneh . The 
SoNIA layouts are created by anchoring the initial node positions, as described 
in Section 14.21 To provide a fair comparison with DMDS and SoNIA, which are 
on-line methods, we use the on-line modification of the Visone method described 



in Sec tion 14.21 DGLL is compared to the CC DR method of ICosta and Hero III 
(120051) . the dynamic spectral layout method of iBrandes et al. I (l2006h (referred to 
here as BFP), and the standard spectral GLL solution dKoren . 2003% 

Summary statistics from the experiments are presented in Tables[T]and|2]for the 
MDS- and GLL-based methods, respectively, and are discussed in the individual 
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Experiment 


Algorithm 


MDS stress 


Centroid cost 


Temporal cost 


Iterations 




DMDS (known) 


0.160 ±0.000 


0.257 ±0.001 


0.262 ±0.001 


45.6 ±0.3 


SBM 


DMDS (learned) 


0.160 ±0.000 


0.308 ± 0.003 


0.293 ± 0.002 


46.9 ± 0.4 


Visone 


0.157 ±0.000 


0.434 ± 0.003 


0.340 ± 0.002 


51.0 ±0.4 




SoNlA 


0.132 ± 0.000 


0.623 ± 0.004 


1.271 ±0.010 


112.9 ± 1.0 




DMDS (learned) 


0.136 


0.656 


0.089 


13.8 


Newcomb 


Visone 


0.107 


1.275 


0.125 


16.9 




SoNIA 


0.065 


1.611 


1.368 


68.5 




DMDS (known) 


0.154 


1.334 


0.207 


37.7 


MIT 


DMDS (learned) 


0.154 


1.491 


0.241 


38.6 


Visone 


0.142 


1.900 


0.290 


45.8 




SoNlA 


0.092 


2.637 


3.384 


108.3 



Table 1: Mean costs of MDS-based layouts (± standard error for SBM simulation experiment). The smallest quantity (within 
one standard error) in each column for each experiment is bolded. DMDS results using both a priori known groups (when 
available) and groups learned by clustering are shown. 
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BFP 


0.636 ± 0.003 


0.637 ± 0.006 


2.269 ± 0.027 




Spectral 


0.605 ± 0.003 


0.945 ± 0.007 


3.179 ±0.022 




DGLL (learned) 


0.820 


1.325 


0.379 


Newcomb 


CCDR 


0.786 


1.334 


1.352 




BFP 


0.783 


1.321 


0.793 




Spectral 


0.761 


1.373 


1.425 




DGLL (known) 


0.130 


1.228 


0.263 




DGLL (learned) 


0.128 


1.357 


0.313 


MIT 


CCDR 


0.099 


1.300 


1.692 




BFP 


0.104 


1.471 


1.897 




Spectral 


0.090 


1.660 


2.478 



Table 2: Mean costs of GLL-based layouts (it standard error for SBM simulation experiment). The smallest quantity (within 
one standard error) in each column for each experiment is bolded. DGLL results using both a priori known groups (when 
available) and groups learned by clustering are shown. 



subsections. The KK choice of MDS weights is used for all of the MDS-based 
methods, and degree-normalized layout is used for all of the GLL-based methods. 

We define three measures of layout quality: static cost, centroid cost, and tem- 
poral cost. The static cost measures how well the current layout coordinates fit the 
current graph snapshot. It is the cost function that would be optimized by the static 
graph layout algorithm, either MDS or GLL. The static cost for the MDS-based 
methods is taken to be the static MDS stress defined in ([T]). The static cost for the 
GLL-based methods is the GLL energy defined in 

The centroid cost is the sum of squared distances between each node and its 
group centroid, which is also the cost function of the well-known method of k- 
means clustering. It is used to measure how close nodes are to members of their 
groupH. When prior knowledge of the groups is available, we calculate the cen- 
troid cost with respect to the known groups, even for the layouts where groups are 
learned by clustering. When prior knowledge is not available, we calculate the 
centroid cost with respect to the learned groups. 

The temporal cost ([T4l ) is the sum of squared distances between node positions 
in layouts at consecutive time ste ps. It is oft en used to quantify how well the 



ment a l map is preserved over t ime ( Brandes and Wagneiill997l : lBaur and SchankL 



20081 : iFrishman and Tall |2008|). As mentioned in Section 13. 1[ the temporal cost 
should only be interpreted as a measure of stability and not a measure of temporal 
goodness-of-fit. 

The costs displayed are appropriately normalized (either by the number of 
nodes or pairs of nodes, depending on the quantity) so they are comparable across 
different data sets. For the MDS-based methods, we also compare the number 
of iterations required for convergence to a tolerance of e = 10^^. For the BFP 
method, the parameter A lies on a different scale from the parameters a and /3 
for the other methods. To ensure a fair comparison, we choose A to minimize 
Cstatic + "Ccentroid + /3Ctemporai> whcrc a and /3 are chosen to be the same parameters 
used for the other methods. 

From Tables [T] and |2l one can see that DMDS and DGLL have lower cen- 
troid and temporal costs than the competing methods in all but one instance due to 
the incorporation of both grouping and temporal regularization. The BFP method 
slightly outperforms DGLL in terms of centroid cost in the Newcomb experiment; 
however, BFP performs comparatively poorly in temporal cost. As expected, the 
lower centroid and temporal costs for DMDS and DGLL are achieved by choosing 
node positions with a higher static cost. Notice also that DMDS requires signif- 
icantly less iterations to converge than SoNIA, which employs no regularization 



^Note that we cannot simply use the grouping cost ( |13t because it is not defined for methods that 
do not incorporate grouping regularization. 
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at all, and slightly less than Visone, which employs only temporal regularization. 
This is an added benefit of regularization. The results for each experiment will be 
discussed in greater detail in the following. 



5.1 Stochastic block model 

In this e xperiment, we genera te simulated networks using a stochastic block model 



(SBM) ([Holland et al.l 119831) . An SBM creates networks with k groups, where 
nodes in a group are stochastically equivalent, i.e. the probability of forming an 
edge between nodes i and j is dependent only on the groups to which i and j 
belong. An SBM is completely specified by the set of probabilities {pcd] c = 
1,...,A:; d = c,c + 1, . . . ,k}, which represent the probability of forming an edge 
between any particular node in group c and any particular node in group d. 

We generate 20 independent samples from a 30-node 4-group SBM with pa- 
rameters Pa = 0.6 and pij = 0.2 for alH / j. Each sample corresponds to a graph 
snapshot at a single time step. The group memberships are randomly assigned at 
the initial time step and remain unchanged up to t = 9. Kit = 10, 1/4 of the nodes 
are randomly re-assigned to different groups to simulate a change in the network 
structure. The group memberships are then held constant until the last time step. 
We create layouts of the network using parameters a = /3 = 1. 

In Fig. m we plot the variation over time of the static, centroid, and tempo- 
ral costs of the MDS-based methods. The costs are averaged over 100 simulation 
runs. As expected, the static cost is higher for the regularized layouts than for 
the unregularized SoNIA layout. The grouping regularization in DMDS results in 
lower centroid cost as expected. When the groups are learned by clustering, the 
centroid cost is slightly higher than with the known groups, but still much lower 
than that of Visone and SoNIA. Although Visone has only temporal regularization, 
notice that it also has a lower centroid cost than SoNIA. This is because the SBM 
parameters are held constant from time steps to 9 and from time steps 10 to 19, so 
that the group structure can be partially revealed by temporal regularization alone 
once enough time samples have been collected. The temporal regularization of 
both DMDS and Visone results in a significantly lower temporal cost than SoNIA. 
The grouping regularization of DMDS also decreases the temporal cost slightly 
compared to Visone. An added benefit of the regularization in DMDS is the sig- 
nificant reduction in the number of iterations required for the MDS algorithm to 
converge, as shown in Table [T] On average, DMDS required less than half as many 
iterations as SoNIA, and slightly less than Visone. 

In Fig. m we plot the variation over time of the static, centroid, and temporal 
costs of the GLL-based methods. Similar to the MDS-based methods, the static 
cost is higher for the regularized layouts, but the centroid and temporal costs are 
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Figure 3: Costs of MDS-based layouts in the SBM experiment at each time step. 
The DMDS layouts have the lowest centroid and temporal costs but also the highest 
MDS stress due to the regularization. 



much lower. Notice that only DGLL is able to generate layouts with low temporal 
cost. The grouping regularization in CCDR reduces the centroid cost but only 
slightly improves the temporal cost. The BFP method, which combines Laplacian 
matrices from two time steps, performs better than the standard spectral method 
both in centroid and temporal cost, but is worse than DGLL in both. 

Notice from Figs. [3] and |4] that the centroid and temporal costs of DMDS and 
DGLL increase at t = 10, reflecting the presence of the change in network struc- 
ture. Such an increase is beneficial for two reasons. First, it indicates that the 
grouping and temporal penalties are not so strong that they mask changes in net- 
work structure, which would be undesirable. Second, it suggests that the centroid 
and temporal costs can be used to filter the sequence of networks for important 
events, such as change points, which can be useful for exploratory analysis of dy- 
namic networks over long periods of time. 

We demonstrate the effect of varying the regularization parameters in DMDS 
and DGLL, respectively, in Figs. [5] and |6] We generate layouts using 10 choices 
each of a and /?, uniformly distributed on a log scale between 0.1 and 10. The 
observations are similar for both DMDS and DGLL. As expected, the temporal 
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Figure 4: Costs of GLL-based layouts in the SBM experiment at each time step. 
The DGLL layouts have the lowest centroid and temporal costs but also the highest 
GLL energy due to the regularization. 



cost decreases for increasing /3. For low values of /3, increasing a also decreases 
the temporal cost. This is a sensible result because nodes can move significantly 
over time but must remain close to the group representative, which lowers the tem- 
poral cost. The result is slightly different when it comes to the centroid cost. As 
expected, increasing a decreases centroid cost. For low values of a, increasing /3 
also decreases centroid cost to a point, but a very high /3 may actually increase cen- 
troid cost, especially in DMDS. This is also a sensible result because a very high 
/3 places too much weight on the initial time step and prevents nodes from moving 
towards their group representative at future time steps. 

From this experiment we can see that there is a coupled effect between group- 
ing and temporal regularization, and that a combination of both can sometimes 
result in better performance with respect to both centroid and temporal cost. How- 
ever, it is important to note that this is not always true. For example, if a node 
changes group between two time steps, then the two penalties can oppose each 
other, with the temporal penalty attempting to pull the node towards its previous 
position and the grouping penalty attempting to pull the node towards its current 
representative, which could be quite far from the node's previous position. This is 
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Figure 5 : Mean centroid and temporal costs of DMDS layouts in the SBM experi- 
ment as functions of a and /3. The centroid and temporal costs decrease as a and 
/3 are increased, respectively; however, a also affects the temporal cost, and /? also 
affects the centroid cost. 




Figure 6: Mean centroid and temporal costs of DGLL layouts in the SBM experi- 
ment as functions of a and /3. The behavior of both costs as functions of a and /3 
is similar to their behavior in DMDS. 



another reason for the increase in both centroid and temporal costs at t = 10, when 
the group structure is altered, in Figs. [3] and ID 



5.2 Newcomb's fraternity 



This data set was collected by Nordlie and Newcomb ((NordlieL 1 19581 : iNewcomb . 
196 ih as part of an experiment on interpersonal relations. It has been examined in 
sever al previous studies including (|Moody et al.Ll2005l : lBender-deMoU and McFarland . 
200a). 17 incoming male transfer students at the University of Michigan were 
housed together in fraternity housing. Each week, the participants ranked their 
preference of each of the other individuals in the house, in private, from 1 to 16. 
Data was collected over 15 weeks in a semester, with one week of data missing 
during week 9, corresponding to Fall break. 



We process the rank data in the same manner as iMoody et all (|2005h and 
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Figure 7: Time plots of 1-D DGLL layouts of Newcomb's fraternity, colored by 
learned groups. Node positions in the layout are relatively stable over time unless 
nodes are changing group. 



Bender-deMoU and McFarland ( 20061) to give a fair comparison with SoNIA. Graph 



snapshots are created by connecting each participant to his 4 most preferred stu- 
dents with weights from 4 decreasing to 1 corresponding to the most preferred 
to the 4th most preferred student. The graph is converted to an undirected graph 
by taking the edge weight between i and j to be the larger of the directed edge 
weights. The weights are converted into dissimilarities for the MDS -based meth- 
ods by dividing each similarity weight by the maximum similarity of 4. No group 
information is known a priori, so the group structure is learned using the AFFECT 
clustering algorithm. 

In Fig. Ul we show a time plot of 1-D layouts created using DGLL, where the 
color of a line segment between time steps t and t + 1 denotes the group mem- 
bership of the node at time step t, and the location of the endpoints correspond 
to the node's position in the layouts at time steps t and t + I. The regularization 
parameters are chosen to be a = /3 = 1. While a 1-D layout does a poor job of 
conveying the topology of the network, some temporal trends can be seen. For 
example, two mostly stable groups form after several weeks, but two nodes appear 
to switch from the red to the blue group around Fall break. 

In Figs. [SlITOl we present a comparison of the first four snapshots from the 
layouts created using DMDS, Visone, and SoNIA, respectively. In both figures, 
the top row corresponds to the layouts, and the bottom row illustrates the node 
movements of each node over time. In the plots on the bottom row, each node is 
drawn twice: once at its current position at time t and once at its previous position 
at time t — 1. An edge connects these two positions; the length of the edge indicates 
how far a node has moved between time steps t — 1 and t. 

At t = 0, the red and green groups are mixed together in the Visone and So- 
NIA layouts, while they are easily distinguished in the DMDS layout due to the 
grouping regularization. Furthermore, the node movements over time in the So- 
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Figure 8: Layouts of Newcomb's fraternity at four time steps (top row) generated 
using proposed DMDS algorithm and node movements between layouts (bottom 
row). The groups remain well-separated. 
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Figure 9: Layouts of Newcomb's fraternity at four time steps (top row) using the 
Visone algorithm and node movements between layouts (bottom row). The groups 
are not as well-separated as in the DMDS layouts. 
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Figure 10: Layouts of Newcomb's fraternity at four time steps (top row) using 
the SoNIA algorithm and node movements between layouts (bottom row). There 
is excessive node movement, and the groups are not as well-separated as in the 
DMDS layouts. 



NIA layouts are much more extreme than in the DMDS layouts. The excessive 
node movement is even more vi sible by compar ing the DMDS and SoNIA anima- 



tions on the supporting website (IXu et al.U2012l) . The excessive node movement is 



reflected in the substantially higher temporal cost of SoNIA compared to DMDS, 
as shown in Table [T] DMDS also outperforms Visone in both mean centroid and 
temporal cost, although the improvement in temporal cost is smaller than compared 
to SoNIA because Visone also employs temporal regularization. Finally, the mean 
number of iterations required for the SoNIA layout to converge is almost four times 
that of DMDS, so DMDS presents significant computational savings in addition to 
better preservation of the mental map. 

For the GL L-based layouts, which can be found on the supporting website 
( Xu et al.[ 2012 ). there is not much difference in terms of the centroid cost. Notice 
from Table |2l that the mean centroid cost differs by only 4% between the best 
(BFP) and worst (spectral). This is not surprising, as the groups are learned using 
evolutionary spectral clustering, which is closely related to GLL. However, DGLL 
performs significantly better than the other methods in terms of temporal cost. The 
temporal smoothing of the Laplacian for BFP helps to lower the temporal cost 
slightly compared to the layouts without any temporal regularization, but it is not 
as effective as the temporal regularization in DGLL. 
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Figure 1 1 : DMDS layouts of MIT Reality Mining data at four time steps using 
the known groups (top row) and node movements between layouts (bottom row). 
Blue nodes denote colleagues working in the same building, and red nodes denote 
incoming students. The incoming students separate from the others after the first 
week of classes (t = 5). 



5.3 MIT Reality Mining 

The MIT Reality Mining data set feagle et aL . 20091) was collected as part of an 
experiment on inferring social networks by using cell phones as sensors. 94 stu- 
dents and staff at MIT were given access to smart phones that were monitored 
over two semesters. The phones were equipped with Bluetooth sensors, and each 
phone recorded the Media Access Control addresses of nearby Bluetooth devices 
at five-minute intervals. Using this proximity data, we construct a sequence of 
graph snapshots where each participant is connected to the 5 participants he or 
she was in highest proximity to during a time step. We divide the data into time 
steps of one w eek, resulting in 46 time steps between August 2004 and June 2005. 
From the MIT academic calendarl . we know the dates of important events such as 
the beginning and end of school terms. We also know that 26 of the participants 
were incoming students at the university's business school, while the rest were col- 
leagues working in the same building. These affiliations are used as the known 
groups. 

The DMDS layouts at three time steps computed using the using the known 
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t = 6 



t = 7 




Figure 12: DMDS layouts of MIT Reality Mining data at four time steps with 
a = 5, /? = 3 using groups learned by clustering (top row) and node movements 
between layouts (bottom row). Colors correspond to learned groups. There is a 
lot of node movement between groups but very little movement within groups, 
resulting in high MDS stress. 



groups with a = 1, /? = 3 are shown in Fig. [TT] A higher value of /3 is chosen 
compared to the previous experiments in order to create more stable layouts due 
to the higher number of nodes. Node labels are not displayed to reduce clutter in 
the figure . We e ncourage readers to view the animation on the supporting website 



(IXu et al.l 12012 ) to get a better idea of the temporal evolution of this network. 



t = 5 corresponds to the first week of classes. Notice that the two groups are 
slightly overlapped at this time step. As time progresses, the group of incoming 
students separates quite clearly from the colleagues working in the same building. 
This result suggests that the incoming students are spending more time in proximity 
with each other than with the remaining participants, which one would expect as 
the students gain familiarity with each other as the semester unfolds. 

The same observation can be made from the DMDS layouts computed using 
the groups learned by the AFFECT clustering algorithm. Initially at t = 5, the 
separation between groups is not clear so many nodes are not correctly classified, 
but at subsequent time steps when the separation is clearer, almost all of the nodes 
are correctly classified. Between time steps 5 and 6 many nodes switch groups. 
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Figure 13: DMDS layouts of MIT Reality Mining data at four time steps with 
a = 1/5, /3 = 3 using groups learned by clustering (top row) and node movements 
between layouts (bottom row). Colors correspond to learned groups. There is more 
movement within groups, resulting in lower MDS stress, but it is more difficult to 
identify movement between groups. 

This can be seen in Fig. [I2j where the colors correspond to the learned groups 
rather than the known groups. These layouts are created using a = 5, /3 = 3; the 
high value of a emphasizes the node movements between groups while sacrificing 
the quality of movements within groups, as discussed in Section 13.41 Notice that 
the groups are very compact and well-separated, so that nodes switching from one 
group to another have large movements between layouts. Compare these layouts 
to those shown in Fig. [131 which are created using a = 1/5, /3 = 3. The lower 
value of a better shows movements within groups, but the large changes in groups 
between time steps 5 and 6 ai^e not as obvious as in Fig. [12] Both layouts are useful 
and provide different insights into the network dynamics; however, the observation 
of the incoming students separating from the other participants is evident in both 
visualizations. 

The benefits of the regularization can be seen once again from the statistics in 
Tables [T] and [2] With group information provided, the DMDS and DGLL layouts 
have lower mean centroid and temporal costs compared to all competing layouts. 
Without group information, the DMDS and DGLL layouts using groups learned by 
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clustering still outperform competing methods in temporal cost, and only CCDR, 
which uses the known group information, outperforms DGLL without group infor- 
mation in terms of centroid cost. The DMDS methods also converge more quickly 
than both Visone and SoNIA. 



6 Conclusions 

In this paper we proposed a regularized graph layout framework for dynamic net- 
work visualization. The proposed framework incorporates grouping and temporal 
regularization into graph layout in order to discourage nodes from deviating too far 
from other nodes in the same group and from their previous position, respectively. 
The layouts are generated in an on-line manner using only present and past data. 
We introduced two dynamic layout algorithms, DMDS and DGLL, which are reg- 
ularized versions of their static counterparts, and demonstrated the importance of 
the regularizers for preserving the mental map in multiple experiments. The pro- 
posed methods generalize existing approaches for temporal regularization in MDS 
and grouping regularization in GLL. 

An important area for future work concerns visualization of extremely large 
dynamic networks containing upwards of thousands of nodes. One issue is related 
to scalability; both DMDS and DGLL require O(n^) computational time and thus 
may not be applicable to extremely large networks. An even more significant issue 
involves the interpretation of visualizations of such large networks. Even when 
equipped with temporal and grouping regularization, layouts of extremely large 
networks may be confusing for a human to interpret so additional techniques may 
be necessary to deal with this challenge. 



A DGLL solution in 2-D 



We derive the expressions for V/, g, H, and J in 2-D. These vectors and ma- 
trices are computed at each iteration in the D GLL algorithm to solve (1251 ) using 
the interior-point algorithm (IByrd et all. Il999h as discussed in Section 13.31 The 
constraints can be written as g{X) = where 
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The gradient of the objective function is given by 

'{2L + 2l3E)9ii - 2f3E9ii[t 



V/(X) 



(2L + 2/3^)^2 - 2/3Ei2[t 
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The Jacobian of the constraints is given by 



J{X) 



'2±lM 




2xfM 



Finally, the Hessian is obtained by 

'2L + 2f5E + 2^iM fisM 



2L + 2/3£; + 2fi2M 



B AFFECT evolutionary clustering algorithm 

Evolutionary clustering algorithms are designed to cluster dynamic data where a 
set of objects is observed over multiple time steps. In the dynamic network set- 
ting, objects correspond to nodes, and observations correspond to graph adjacency 
matrices W[t]. The AFFECT (A daptive Forgetting Factor for Evolutionary Clus- 
tering and Tracking) framework dXu et all 12011 a ) involves creating a smoothed 
adjacency matrix at each time step and then performing ordinary static clustering 
on this matrix. The smoothed adjacency matrix is given by 

■^[t] = a[t]^[t - 1] + (1 - a[t\)W[t\, 

where a\t] is a forgetting factor that controls how quickly previous adjacency ma- 
trices are forgotten. We drop the time index for quantities at time t for simplicity. 
The AFFECT framework adaptively estimates the optimal amount of smoothing to 
apply at each time step in order to minimize mean-squared error (MSE) in terms of 
the Frobenius norm E[||^ — ^'|||'], where ^ denotes the expected adjacency ma- 
trix E[M^], which can be viewed as a ma trix of unknown states characterizing the 
network structure at time t. It is shown in lxu et al. ( 2011a ) that the optimal choice 
of a is given by 



a 



En \-^n [ 
i=iE,=iVar(m 



YTi=\ YTj=\ \ - 1] - V'ii ) + var (wij, 



For real networks, and var(ii;jj) are unknown so a* cannot be computed. 

The AFFECT framework iteratively estimates the optimal forgetting factor and 
clusters nodes using a two-step procedure. Begin by initializing the clusters; for 
example, using the clusters at the previous time step, a* is estimated by replacing 
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the unknown means and variances with sample means and variances computed over 
the clusters. Static clustering is then performed on the smoothed adjacency matrix 
^' to obtain cluster memberships. The estimate of a* and the cluster memberships 
are then refined by iterating t he tw o steps. In this paper, we use normalized cut 
spectral clustering ( Ng et all 2001 ) as the static clustering algorithm. We refer 
interested readers to lXu et al.l (120 Hal ) for additional details. 
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